import numpy as np


# Gauss列主消元法
# 输入：矩阵a及其行/列数n, 矩阵b
def Gauss(n, a, b):
    m = np.zeros((n, n))    # 生成一个大小为n*n的全零矩阵
    x = np.zeros(n)         # 生成一个大小为n的全零向量

   # 消元过程
    for k in range(n - 1):
        max, p = abs(a[k][k]), k

        # 寻找主元并标记主元位置
        for j in range(k + 1, n):
            if(abs(a[j][k]) > max):
                max = abs(a[j][k])
                p = j

        # 主元为0时返回奇异
        if(a[p][k] == 0):
            print("Singular!")
            return

        # 交换当前行和主元所在行
        elif(p != k):
            for j in range(n):
                a[p][j], a[k][j] = a[k][j], a[p][j]
            b[k], b[p] = b[p], b[k]

        # 消元
        for i in range(k + 1, n):
            m[i][k] = a[i][k] / a[k][k]
            for j in range(n):
                a[i][j] = a[i][j] - a[k][j] * m[i][k]
            b[i] = b[i] - b[k] * m[i][k]

    # 回代过程
    if(a[n - 1][n - 1] == 0):
        print("Singular!")
        return
    x[n - 1] = float(format(b[n - 1] / a[n - 1][n - 1], '.8f'))
    for k in range(n - 2, -1, -1):
        tot = 0
        for j in range(k + 1, n):
            tot = tot + a[k][j] * x[j]
        x[k] = float(format((b[k] - tot) / a[k][k], '.8f'))

    # 打印结果
    print("x = ")
    for k in range(n):
        print(x[k])
    return


# 打印题号
print("-----T1(1)-----")
a = np.array([[0.4096, 0.1234, 0.3678, 0.2943], [0.2246, 0.3872, 0.4015, 0.1129], [0.3645, 0.1920, 0.3781, 0.0643],[0.1784, 0.4002, 0.2786, 0.3927]])
b = [1.1951, 1.1262, 0.9989, 1.2499]
Gauss(4, a, b)


print("-----T1(2)-----")
a = np.array([[136.01, 90.860, 0.0, 0.0], [90.860, 98.810, -67.590, 0.0], [0.0, -67.590, 132.01, 46.260], [0.0, 0.0, 46.260, 177.17]])
b = [226.87, 122.08, 110.68, 223.43]
Gauss(4, a, b)

print("-----T1(3)-----")
a = np.array([[1, 1/2, 1/3, 1/4], [1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7]])
b = [25/12, 77/60, 57/60, 319/420]
Gauss(4, a, b)

print("-----T1(4)-----")
a = np.array([[10.0, 7.0, 8.0, 7.0], [7.0, 5.0, 6.0, 5.0], [8.0, 6.0, 10.0, 9.0], [7.0, 5.0, 9.0, 10.0]])
b = [32.0, 23.0, 33.0, 31.0]
Gauss(4, a, b)

print("-----T2(1)-----")
a = np.array([[197.0, 305.0, -206.0, -804.0], [46.8, 71.3, -47.4, 52.0], [88.6, 76.4, -10.8, 802.0], [1.45, 5.90, 6.13, 36.5]])
b = [136.0, 11.7, 25.1, 6.60]
Gauss(4, a, b)

print("-----T2(2)-----")
a = np.array([[0.5398, 0.7161, -0.5554, -0.2982], [0.5257, 0.6924, 0.3565, -0.6225], [0.6465, -0.8187, -0.1872, 0.1291], [0.5814, 0.9400, -0.7779, -0.4042]])
b = [0.2058, -0.0503, 0.1070, 0.1859]
Gauss(4, a, b) 

print("-----T2(3)-----")
a = np.array([[10.0, 1.0, 2.0], [1.0, 10.0, 2.0], [1.0, 1.0, 5.0]])
b = [13.0, 13.0, 7.0]
Gauss(3, a, b) 

print("-----T2(4)-----")
a = np.array([[4.0, -2.0, -4.0], [-2.0, 17.0, 10.0], [-4.0, 10.0, 9.0]])
b = [-2.0, 25.0, 15.0]
Gauss(3, a, b) 
